﻿using System;
using SpiceNet.Circuits;
using SpiceNet.Diagnostics;
using SpiceNet.Models;
using System.Numerics;

namespace SpiceNet.Simulations
{
    /// <summary>
    /// A class capable of performing a noise analysis
    /// </summary>
    public class Noise : Analysis
    {
        // Enumerations
        public enum OutputMode { N_DENS = 1, INT_NOIZ = 2 }
        public enum Operation { N_OPEN = 1, N_CALC = 2, N_CLOSE = 3 }
        public enum Sources { SHOTNOISE = 1, THERMNOISE = 2, N_GAIN = 3 }
        public enum StepTypes { DECADE = 1, OCTAVE = 2, LINEAR = 3 }

        // Tolerances
        public const double N_MINLOG = 1e-38;
        public const double N_MINGAIN = 1e-20;
        public const double N_INTFTHRESH = 1e-10;
        public const double N_INTUSELOG = 1e-10;

        // Maximum length of noise output variables we will generate
        public const int N_MXVLNTH = 64;

        /// <summary>
        /// A structure for saving and retrieving integrated noise data
        /// </summary>
        public class InstanceNoiseData
        {
            public double LNLSTDENS = 0.0;
            public double OUTNOIZ = 0.0;
            public double INNOIZ = 0.0;
        }

        /// <summary>
        /// Constructor
        /// </summary>
        public Noise() : base(typeof(NoiseJob))
        {
            Name = "Noise";
            Description = "Noise analysis";
            Domain = DomainTypes.FrequencyDomain;
            DoIc = true;
        }

        /// <summary>
        /// Wrapper for NOISEan in noisean.c
        /// Perform noise analysis
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="reset"></param>
        public override void Func(Circuit ckt, bool reset)
        {
            NoiseJob job = ckt.CurrentJob as NoiseJob;

            Complex Val;
            int posOutNode = 0;
            int negOutNode = 0;
            int step;
            NoiseData data = new NoiseData();
            double freqTol;

            if (ckt.Nodes.ContainsNode(job.output))
                posOutNode = ckt.Nodes[job.output].Number;
            if (ckt.Nodes.ContainsNode(job.outputRef))
                negOutNode = ckt.Nodes[job.outputRef].Number;

            // See if the source specified is AC
            // Find the type codes for voltage sources and current sources
            Model[] vsrcs = ckt.FindModels(typeof(VoltageSourceModel));
            Model[] isrcs = ckt.FindModels(typeof(CurrentSourceModel));

            // See if the source specified is AC
            Instance inst = null;

            // Loop through voltage sources if any
            for (int i = 0; i < (vsrcs?.Length ?? 0); i++)
            {
                inst = vsrcs[i].Find(job.input);
                if (!(inst as VoltageSource)?.VSRCacMag.Given ?? false)
                    throw new CircuitException("Noise input source has no AC value");
            }

            // Loop through current sources if any
            for (int i = 0; i < (isrcs?.Length ?? 0); i++)
            {
                inst = isrcs[i].Find(job.input);
                if (!(inst as CurrentSource)?.ISRCacMag.Given ?? false)
                    throw new CircuitException("Noise input source has no AC value");
            }

            if ((job.NsavFstp == 0) || reset)
            {
                switch (job.NstpType)
                {
                    case StepTypes.DECADE:
                        job.NfreqDelta = Math.Exp(Math.Log(10.0) / job.NnumSteps);
                        break;

                    case StepTypes.OCTAVE:
                        job.NfreqDelta = Math.Exp(Math.Log(2.0) / job.NnumSteps);
                        break;

                    case StepTypes.LINEAR:
                        job.NfreqDelta = (job.NstopFreq - job.NstartFreq) / (job.NnumSteps + 1);
                        break;
                }

                // Calculate the operating point
                Op(ckt,
                    (ckt.Mode & Circuit.Modes.Uic) | Circuit.Modes.DcOp | Circuit.Modes.InitJct,
                    (ckt.Mode & Circuit.Modes.Uic) | Circuit.Modes.DcOp | Circuit.Modes.InitFloat,
                    ckt.Config.DcMaxIterations);

                step = 0;
                data.freq = job.NstartFreq;
                data.outNoiz = 0.0;
                data.inNoise = 0.0;

                // The current front-end needs the namelist to be fully declared before an OUTpBeginplot
                // Not needed here...

                data.numPlots = 0; // We don't have any plots yet
                NoiseLoad(ckt, OutputMode.N_DENS, Operation.N_OPEN, data);

                // All names in the namelist have been declared. Now start the plot
                // Not needed here
            }
            else // We must have paused before. Pick up where we left off
            {
                step = job.NsavFstp;
                switch (job.NstpType)
                {
                    case StepTypes.DECADE:
                    case StepTypes.OCTAVE:
                        data.freq = job.NstartFreq * Math.Exp(step * Math.Log(job.NfreqDelta));
                        break;
                    case StepTypes.LINEAR:
                        data.freq = job.NstartFreq + step * job.NfreqDelta;
                        break;
                    default:
                        throw new CircuitException("Invalid step type");
                }
                job.NsavFstp = 0;
                data.outNoiz = job.NsavOnoise;
                data.inNoise = job.NsavInoise;
            }

            switch (job.NstpType)
            {
                case StepTypes.DECADE:
                case StepTypes.OCTAVE:
                    freqTol = job.NfreqDelta * job.NstopFreq * ckt.Config.RelTol;
                    break;
                case StepTypes.LINEAR:
                    freqTol = job.NfreqDelta * ckt.Config.RelTol;
                    break;
                default:
                    throw new CircuitException("Invalid step type");
            }

            data.lstFreq = data.freq;

            // Do the noise analysis over all frequencies
            while (data.freq <= job.NstopFreq + freqTol)
            {
                // Test for pausing - not implemented yet

                ckt.StateComplex.Omega = new Complex(0.0, 2.0 * Math.PI * data.freq);
                ckt.Mode = (ckt.Mode & Circuit.Modes.Uic) | Circuit.Modes.Ac;

                // Solve the original AC system to get the transfer function between the input and output

                AcIterate(ckt);
                Val = ckt.StateComplex.Solution[posOutNode] - ckt.StateComplex.Solution[negOutNode];
                data.GainSqInv = 1.0 / Math.Max((Val.Real * Val.Real) + (Val.Imaginary * Val.Imaginary), N_MINGAIN);
                data.lnGainInv = Math.Log(data.GainSqInv);

                // Set up a block of "common" data so we don't have to recalculate it for every device
                data.delFreq = data.freq - data.lstFreq;
                data.lnFreq = Math.Log(Math.Max(data.freq, N_MINLOG));
                data.lnLastFreq = Math.Log(Math.Max(data.lstFreq, N_MINLOG));
                data.delLnFreq = data.lnFreq - data.lnLastFreq;

                if ((job.NStpsSm != 0) && ((step % (job.NStpsSm)) == 0))
                    data.prtSummary = true;
                else
                    data.prtSummary = false;

                /*
                data.outNumber = 1;
                 */

                data.outNumber = 0;
                /* The frequency will NOT be stored in array[0] as before; instead,
                 * it will be given in refVal.rValue (see later)
                 */
                NinzIterate(ckt, posOutNode, negOutNode); // Solve the adjoint system

                /* Now we use the adjoint system to calculate the noise
                 * contributions of each generator in the circuit
                 */
                NoiseLoad(ckt, OutputMode.N_DENS, Operation.N_CALC, data);
                if (data.prtSummary)
                    Dump(ckt);
                data.lstFreq = data.freq;

                // Update the frequency
                switch (job.NstpType)
                {
                    case StepTypes.DECADE:
                    case StepTypes.OCTAVE:
                        data.freq *= job.NfreqDelta;
                        break;
                    case StepTypes.LINEAR:
                        data.freq += job.NfreqDelta;
                        break;
                    default:
                        throw new CircuitException("Wrong step type");
                }
                step++;
            }

            NoiseLoad(ckt, OutputMode.N_DENS, Operation.N_CLOSE, data);

            data.numPlots = 0;
            data.outNumber = 0;

            if (job.NstartFreq != job.NstopFreq)
            {
                NoiseLoad(ckt, OutputMode.INT_NOIZ, Operation.N_OPEN, data);

                // Begin a plot for the integrated noise
                NoiseLoad(ckt, OutputMode.INT_NOIZ, Operation.N_CALC, data);

                // Do exporting data
                Dump(ckt);

                NoiseLoad(ckt, OutputMode.INT_NOIZ, Operation.N_CLOSE, data);
            }
        }

        /// <summary>
        /// Load AC version
        /// </summary>
        public void AcLoad(Circuit ckt)
        {
            // Clear the state
            ckt.StateComplex.Matrix.Clear();
            ckt.StateComplex.Rhs.Clear();

            // Iterate through all models
            foreach (Type t in ckt.ModelTypes)
            {
                var models = ckt[t];
                for (int i = 0; i < models.Length; i++)
                    models[i].AcLoad(ckt);
            }
        }

        /// <summary>
        /// Iterate for complex analysis
        /// </summary>
        public void AcIterate(Circuit ckt)
        {
            // retry:
            ckt.IsCon = true;

            // Load AC
            AcLoad(ckt);

            // Solve
            // OldSolution = Solution;
            ckt.StateComplex.Solve();
        }

        /// <summary>
        /// Iterate for noise analysis
        /// </summary>
        /// <param name="posDrive"></param>
        /// <param name="negDrive"></param>
        public void NinzIterate(Circuit ckt, int posDrive, int negDrive)
        {
            // Clear out the right hand side vector
            ckt.StateComplex.Rhs.Clear();

            // Apply unit current excitation
            ckt.StateComplex.Rhs[posDrive] = 1.0;
            ckt.StateComplex.Rhs[negDrive] = -1.0;

            // Solve (note: index 0 is already set to 0 + 0i by construction)
            // Solve the transposed matrix
            ckt.StateComplex.SolveTransposed();
        }

        /// <summary>
        /// Load the structures with noisedata
        /// </summary>
        public void NoiseLoad(Circuit ckt, OutputMode mode, Operation operation, NoiseData data)
        {
            double outNdens = 0.0;
            ckt.CurrentNoiseData = data;

            // Load the noise data for each instance
            foreach (Type t in ckt.ModelTypes)
            {
                var models = ckt[t];
                for (int i = 0; i < models.Length; i++)
                    models[i].NoiseLoad(ckt, mode, operation, data, ref outNdens);
            }

            // Decide how to fill the NoiseData based on the operation and mode
            switch (operation)
            {
                case Operation.N_OPEN:

                    // Take care of the noise for the circuit as a whole
                    switch (mode)
                    {
                        case OutputMode.N_DENS:

                            // Add onoise_spectrum and inoise_spectrum
                            data.namelist.Add("onoise_spectrum");
                            data.namelist.Add("inoise_spectrum");

                            // Allocate the output vector
                            data.outpVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(data.namelist.Count);
                            break;

                        case OutputMode.INT_NOIZ:

                            // Add onoise_total and inoise_total
                            data.namelist.Add("onoise_total");
                            data.namelist.Add("inoise_total");

                            // Allocate the output vector
                            data.outpVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(data.namelist.Count);
                            break;

                        default:
                            break;
                    }
                    break;

                case Operation.N_CALC:

                    switch (mode)
                    {
                        case OutputMode.N_DENS:
                            if (ckt.CurrentJob == null || (!(ckt.CurrentJob is NoiseJob)))
                                throw new CircuitException("No job specified");
                             
                            if (((ckt.CurrentJob as NoiseJob).NstpType == 0) || (data.prtSummary))
                            {
                                data.outpVector[data.outNumber++] = outNdens;
                                data.outpVector[data.outNumber++] = (outNdens * data.GainSqInv);

                                // We are handling output data a bit differently (see the noise analysis)
                            }
                            break;
                        case OutputMode.INT_NOIZ:
                            data.outpVector[data.outNumber++] = data.outNoiz;
                            data.outpVector[data.outNumber++] = data.inNoise;
                            break;
                        default:
                            break;
                    }
                    break;

                case Operation.N_CLOSE:
                    // Closing the plot, but this is not needed here...
                    break;

            }
        }
    }
}
